1 Load packages

library("knitr")       # for knitting 
library("janitor")     # for cleaning variable names 
library("RSQLite")     # for reading in db files 
library("tidyjson")    # for parsing json files 
library("corrr")       # for correlations
library("DT")          # tables in Rmarkdown
library("grid")        # functions for dealing with images 
library("xtable")      # for latex tables
library("png")         # adding pngs to images
library("egg")         # for geom_custom
library("patchwork")   # for making figure panels
library("ggtext")      # for text in ggplots 
library("kableExtra")  # for knitr tables
library("lubridate")   # for dealing with dates
library("modelr")      # for bootstrapping
library("rsample")     # for bootstrapping
library("tidyverse")   # for everything else 
theme_set(theme_classic() + 
    theme(text = element_text(size = 24)))

opts_chunk$set(comment = "",
               fig.show = "hold")

options(dplyr.summarise.inform = F)

# root mean squared error
rmse = function(x, y){
  return(sqrt(mean((x - y)^2)))
}

# function for printing out html or latex tables 
print_table = function(data, format = "html", digits = 2){
  if(format == "html"){
    data %>% 
      kable(digits = digits) %>% 
      kable_styling()
  }else if(format == "latex"){
    data %>% 
      xtable(digits = digits) %>%
      print(include.rownames = F,
            booktabs = T)
  }
}

2 DATA

con = dbConnect(SQLite(),dbname = "../../data/participants_anonymized.db");
df = dbReadTable(con,"inference_from_blame")
dbDisconnect(con)

df = df %>% 
  filter(status %in% c(4, 5))

2.1 Trial info

# trial information 
df.trialinfo = df$datastring %>% 
  as.tbl_json() %>% 
  enter_object("data") %>% 
  gather_array() %>% 
  enter_object("trialdata") %>% 
  gather_object("name") %>% 
  filter(document.id == 1,
         name == "situation") %>% 
  gather_object("index") %>% 
  append_values_string() %>% 
  as_tibble() %>% 
  filter(index != "unknowns",) %>% 
  pivot_wider(names_from = index, 
              values_from = string) %>% 
  select(trial = id,
         situation = base_image) %>% 
  mutate(situation = str_remove(situation, ".png")) %>% 
  filter(!trial %in% c("att1", "att2")) %>% 
  separate(situation,
           into = c("trees", "a", "b", "c")) %>% 
  mutate(trees = as.numeric(str_remove(trees, "t")),
         trial = as.numeric(trial) + 1,
         across(a:c,
                ~ str_extract(., "\\-*\\d+\\.*\\d*"),
                .names = "{.col}_strength"),
         across(a:c,
                ~ str_extract(., "fish|None|trees"),
                .names = "{.col}_choice"),
         across(a:c,
                ~ str_extract(., "low|medium|high"),
                .names = "{.col}_blame"),
         ) %>% 
  select(-c(a:c)) %>% 
  arrange(trial)

2.2 Model predictions by trial

df.model = read_csv("../../data/model_results/model_posteriors.csv") %>% 
  rename(trial = trial_id) %>% 
  rename_with(.fn = ~ str_c(c("a", "b", "c"),
                            rep(c("_strength", "_choice", "_blame"), each = 3)),
              .cols = contains("_")) %>% 
  mutate(trial = trial + 1,
         across(contains("choice"), ~ ifelse(. == 0, "fish", "trees"))) %>% 
  select(-contains("blame")) %>% 
  pivot_wider(names_from = model,
              values_from = post) %>% 
  relocate(trial) #%>% 
  #select(-reward_maximizer)

df.model.long = df.model %>%
  pivot_longer(cols = -c(trial, rationality, mixture, pivotality, random),
               names_to = "question",
               values_to = "choice",
               values_transform = list(choice = as.character)) %>%
  relocate(question, choice, .after = trial) %>%
  arrange(trial, question, choice)

2.3 Main info

# main data 
df.long = df$datastring %>% 
  as.tbl_json() %>% 
  enter_object("data") %>% 
  gather_array() %>% 
  enter_object("trialdata") %>% 
  gather_object() %>% 
  append_values_string() %>% 
  as_tibble() %>% 
  rename(participant = document.id,
         trial = array.index) %>% 
  mutate(test = string == "TESTTRIAL") %>% 
  group_by(participant, trial) %>% 
  filter(any(test)) %>% 
  ungroup() %>% 
  select(participant, trial, name, string) %>% 
  pivot_wider(names_from = name, 
              values_from = string) %>% 
  clean_names() %>% 
  select(participant, -trial, trial = trial_ix, order = trial_number, a_strength:fish) %>% 
  mutate(across(a_strength:fish, ~ na_if(., "???")))

# filter out participants who didn't pass both of the attention checks 
df.long = df.long %>% 
  mutate(include = trial == "att1" & fish == "4",
         include = ifelse(trial == "att2" & fish == "0", T, include)) %>% 
  group_by(participant) %>% 
  filter(sum(include) == 2) %>% 
  ungroup() %>% 
  filter(!trial %in% c("att1", "att2"))

# some more restructuring 
df.long = df.long %>% 
  mutate(trial = as.numeric(trial) + 1,
         across(.cols = contains("strength"),
                .fns = ~ factor(., levels = 1:3)),
         across(.cols = contains("choice"),
                .fns = ~ factor(., levels = c("fish", "trees"))),
         trees = str_extract(trees, "\\-*\\d+\\.*\\d*"),
         trees = factor(trees, levels = 1:3)) %>% 
  select(-c(fish, order, include)) %>% 
  pivot_longer(cols = -c(participant, trial),
               names_to = "question",
               values_to = "choice") %>% 
  na.omit() %>% 
  arrange(participant, trial)

# writing out for use in python
df.long %>% write_csv('../../data/df.long.csv')

# aggregated results combined with model predictions 
df.aggregate = df.long %>% 
  count(trial, question, choice) %>% 
  group_by(trial, question) %>% 
  mutate(people = n/sum(n),
         choice = as.character(choice)) %>% 
  ungroup() %>% 
  left_join(df.model.long,
            by = c("trial", "question", "choice")) %>% 
  select(-n) %>% 
  # double check that these steps below are valid
  group_by(trial, question, choice) %>% 
  summarize(people = mean(people),
            across(.cols = c(rationality:mixture),
                   .fns = ~ sum(.))) %>% 
  ungroup() %>% 
  pivot_longer(cols = people:mixture,
               names_to = "index",
               values_to = "value")

2.4 Demographic info

# demographics 
df.demographics = df$datastring %>% 
  as.tbl_json() %>%
  enter_object("questiondata") %>% 
  gather_object() %>% 
  append_values_string() %>% 
  as_tibble() %>% 
  rename(participant = document.id) %>% 
  pivot_wider(names_from = name,
              values_from = string) %>% 
  mutate(begin = df$beginhit,
         end =  df$endhit,
         time = as.duration(interval(ymd_hms(df$beginexp), ymd_hms(df$endhit)))) %>% 
  select(-c(begin, end))

df.demographics %>% 
  filter(participant %in% unique(df.long$participant)) %>% 
  mutate(age = as.numeric(age)) %>% 
  summarize(age_mean = mean(age),
            age_sd = sd(age),
            n_female = sum(sex == "female"),
            n_male = sum(sex == "male"),
            n_nosay = sum(sex == "noresponse"),
            time_mean = mean(time) / 60,
            time_sd = sd(time) / 60) %>% 
  print_table(digits = 1)
age_mean age_sd n_female n_male n_nosay time_mean time_sd
40.9 10.8 24 25 1 12 4.1

2.5 Feedback

The question prompt was: “What factors did you take into account when trying to fill in the missing information? Do you have any comments about the experiment?”

datatable(df.demographics %>% 
            select(participant, feedback))

3 STATS

3.1 Bootstrap confidence intervals on choices

set.seed(1)

# percentages with bootstrapped confidence intervals 
df.confidence = df.long %>% 
  mutate(index = str_c(trial, ".", question)) %>%
  group_by(index) %>% 
  nest() %>% 
  mutate(bootstraps = map(.x = data,
                          .f = ~ bootstrap(.x, n = 100))) %>% 
  unnest(bootstraps) %>% 
  mutate(counts = map(.x = strap, 
                      .f = ~ .x %>% 
                        as_tibble() %>% 
                        count(choice))) %>% 
  select(index, .id, counts) %>% 
  unnest(counts) %>% 
  group_by(index, .id) %>% 
  mutate(p = n / sum(n)) %>% 
  group_by(index, choice) %>% 
  summarize(low = quantile(p, probs = 0.025),
            high = quantile(p, probs = 0.975)) %>% 
  ungroup() %>% 
  separate(index,
           into = c("trial", "question"),
           sep = "\\.") %>% 
  mutate(trial = as.numeric(trial))

df.aggregate = df.aggregate %>% 
  left_join(df.confidence,
            by = c("trial", "question", "choice")) %>% 
  mutate(across(c(low, high), ~ ifelse(index != "people", NA, .)))

# save("df.aggregate",
#      file = "../../data/df_aggregate.RData")
load(file = "../../data/df_aggregate.RData")

3.2 Correlation matrix

df.aggregate %>% 
  pivot_wider(names_from = index,
              values_from = value) %>% 
  select(rationality:mixture) %>% 
  correlate() %>% 
  shave() %>% 
  fashion() %>% 
  print_table()
term rationality pivotality mixture
rationality
pivotality .57
mixture .96 .70

3.3 Model summary statistics (best fitting parameters, BIC for each model, n ppt best fit)

# per participant model predictions
df.fit.model = read_csv("../../data/model_results/gen_mix_posts.csv")

df.fit.model =
  df.fit.model %>%
  mutate(model_name = case_when(
    w == 0 ~ "pivotality",
    w == 1 ~ "rationality",
    TRUE ~ "mixture"))

# likelihood of judgments if responding randomly, computed from python
random_loglk = -40.111855503577154

bic = function(k, n, log_lk) (k * log(n)) - (2 * log_lk)

# find best fitting model over all participants
model_summary = df.fit.model %>%
  group_by(w, k, rationality_beta, decision_beta, model_name) %>%
  summarize(loglk = sum(log(lk))) %>%
  mutate(bic_ = case_when(
    model_name == "pivotality" ~ bic(1, 36 * 49, loglk),
    model_name == "rationality" ~ bic(3, 36 * 49, loglk),
    model_name == "mixture" ~ bic(4, 36 * 49, loglk)
  )) %>%
  group_by(model_name) %>%
  mutate(best_fit = min(bic_)) %>%
  filter(best_fit == bic_) %>%
  select(-best_fit) %>%
  select(-loglk) %>%
  bind_rows(tibble(w = NA,
    k = NA,
    rationality_beta = NA,
    decision_beta = NA,
    model_name = "random",
    bic_ = bic(0, 36 * 49, 49 * random_loglk))) %>%
  rename(model = model_name,
         beta_r = rationality_beta,
         beta_d = decision_beta, BIC = bic_) %>%
  relocate(model, .before = w) %>%
  arrange(-BIC)

# getting best fitting by participant
model_summary_ppt = df.fit.model %>%
  group_by(model_name, w, k, rationality_beta, decision_beta, participant) %>%
  summarize(loglk = sum(log(lk))) %>%
  mutate(bic_ = case_when(
    model_name == "pivotality" ~ bic(1, 36, loglk),
    model_name == "rationality" ~ bic(3, 36, loglk),
    model_name == "mixture" ~ bic(4, 36, loglk)
  )) %>%
  ungroup() %>%
  group_by(participant) %>%
  mutate(best_fit = min(bic_)) %>%
  filter(best_fit == bic_) %>%
  select(-best_fit) %>%
  mutate(random = bic(0, 36, random_loglk)) %>%
  mutate(best_model = case_when(
    random < bic_ ~ "random",
    model_name == "pivotality" ~ "pivotality",
    model_name == "rationality" ~ "rationality",
    model_name == "mixture" ~ "mixture"
  )) %>%
  rename(model = best_model)

model_ppt_bestfit = model_summary_ppt %>%
  select(participant, model)

model_summary %>% 
  left_join(model_summary_ppt %>%
              group_by(model) %>%
              tally(),
            by = "model") %>%
  print_table()
model w k beta_r beta_d BIC n
random NA NA NA NA 3930.96 288
pivotality 0.0 1 0.5 1.5 3841.22 90
pivotality 0.0 1 1.0 1.5 3841.22 90
pivotality 0.0 1 1.5 1.5 3841.22 90
pivotality 0.0 1 3.0 1.5 3841.22 90
pivotality 0.0 1 9.0 1.5 3841.22 90
pivotality 0.0 2 0.5 1.5 3841.22 90
pivotality 0.0 2 1.0 1.5 3841.22 90
pivotality 0.0 2 1.5 1.5 3841.22 90
pivotality 0.0 2 3.0 1.5 3841.22 90
pivotality 0.0 2 9.0 1.5 3841.22 90
pivotality 0.0 3 0.5 1.5 3841.22 90
pivotality 0.0 3 1.0 1.5 3841.22 90
pivotality 0.0 3 1.5 1.5 3841.22 90
pivotality 0.0 3 3.0 1.5 3841.22 90
pivotality 0.0 3 9.0 1.5 3841.22 90
rationality 1.0 2 1.0 1.5 3652.89 15
mixture 0.9 2 1.0 1.5 3651.88 7

4 PLOTS

4.1 Scatter plots

df.plot = df.aggregate %>%
  select(-c(low, high)) %>% 
  pivot_wider(names_from = "index",
              values_from = "value") %>% 
  mutate(question_type = str_extract(question, "trees|strength|choice"))

df.models = tibble(name = c("rationality", "pivotality", "mixture"),
                   color = c("blue", "red", "purple"))

fun_scatter = function(df_plot, df_models, index){
  p = ggplot(data = df_plot,
             mapping = aes_string(x = df_models$name[index],
                                  y = "people")) + 
    geom_smooth(method = "lm",
                mapping = aes(group = 1),
                fill = df_models$color[index],
                color = df_models$color[index],
                alpha = 0.1,
                show.legend = F) +  
    geom_point(mapping = aes(fill = question_type),
               size = 2.5,
               shape = 21,
               alpha = 0.75) + 
    annotate(geom = "text",
             label = df_plot %>% 
               summarize(r = cor(people, .data[[df_models$name[index]]]),
                         rmse = rmse(people, .data[[df_models$name[index]]])) %>% 
               mutate(across(.fns = ~ round(., 2))) %>% 
               unlist() %>% 
               str_c(names(.), " = ", .),
             x = c(0, 0), 
             y = c(1, 0.92),
             size = 7,
             hjust = 0) +
    scale_x_continuous(breaks = seq(0, 1, 0.25),
                       limits = c(0, 1)) + 
    scale_y_continuous(breaks = seq(0, 1, 0.25),
                       limits = c(0, 1)) + 
    scale_fill_manual(values = c("gray80", "gray40", "gray0")) + 
    labs(x = df_models$name[index],
         fill = "question") + 
    theme(legend.position = "none",
          axis.title.x = element_markdown(color = df_models$color[index]))
  
  # add legend to second plot
  if (index == 1){
    p = p +
      theme(legend.position = c(1, 0),
            legend.justification = c(1, -0.1)) +
      guides(fill = guide_legend(override.aes = list(size = 6)))
  }
  return(p)
}


l.plots = map(.x = 1:3, ~ fun_scatter(df.plot, df.models, .x))

wrap_plots(l.plots,
           ncol = 3) + 
  plot_annotation(tag_levels = "A")

# ggsave(str_c("../../figures/plots/scatterplots.pdf"),
#        width = 18,
#        height = 6)

4.2 Bar chart and model predictions (all trials)

df.plot = df.aggregate %>% 
  select(-c(low, high)) %>% 
  mutate(choice = str_replace(choice, "trees", "T"),
         choice = str_replace(choice, "fish", "F"),
         item = str_c(question, choice),
         item = str_replace(item, "trees", "t"),
         item = str_remove(item, "_strength"),
         item = str_remove(item, "_choice"),
         question_type = str_extract(question, "trees|strength|choice"))
  
ggplot(data = df.plot,
       mapping = aes(x = item,
                     y = value)) +
  geom_col(data = df.plot %>% 
             filter(index == "people"),
           mapping = aes(fill = question_type,
                         linetype = question_type),
           alpha = 0.5,
           color = "black") +
  geom_point(data = df.plot %>% 
               filter(index %in% c("rationality", "pivotality", 
                                   "mixture")) %>% 
               mutate(index = factor(index, levels = c("rationality", "pivotality", 
                                   "mixture"))),
             mapping = aes(shape = index,
                           fill = index),
             position = position_dodge(width = 0.5),
             color = "black",
             stroke = 1,
             size = 4) +
  facet_wrap(~trial, scales = "free_x") +
  labs(linetype = "question",
       shape = "model") + 
  scale_y_continuous(breaks = seq(0, 1, 0.25),
                     labels = seq(0, 100, 25),
                     expand = expansion(add = c(0, 0.1))) +
  scale_shape_manual(values = c(21, 23, 24)) + 
  scale_linetype_manual(values = rep(1, 3)) + 
  scale_fill_manual(breaks = c("choice", "strength", "trees",
                               "rationality", "pivotality", "mixture"),
                    values = c(choice = "gray80", 
                               strength = "gray40", 
                               trees = "gray0", 
                               rationality = "blue", 
                               pivotality = "red", 
                               mixture = "purple")) + 
  theme(legend.position = "bottom",
        panel.grid.major.y = element_line(),
        axis.text.y = element_text(size = 14),
        axis.title = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        panel.spacing.x = unit(1, "cm"),
        panel.spacing.y = unit(1, "cm")) + 
  guides(fill = F,
         linetype = guide_legend(override.aes = list(fill = c("gray80",
                                                              "gray40",
                                                              "gray0"))),
         shape = guide_legend(override.aes = list(fill = c("blue", "red", "purple"))))

# ggsave("../../figures/plots/all_trials.pdf",
#        width = 20,
#        height = 12)

4.3 Bar chart and model predictions (selection)

# plotting example trials
trials = c(1, 3, 9, 11, 18, 34, 25, 10)
labels = 1:8

df.plot = df.aggregate %>% 
  mutate(choice = str_replace(choice, "trees", "T"),
         choice = str_replace(choice, "fish", "F"),
         item = str_c(question, choice),
         item = str_replace(item, "trees", "t"),
         item = str_remove(item, "_strength"),
         item = str_remove(item, "_choice"),
         question_type = str_extract(question, "trees|strength|choice")) %>%
  filter(trial %in% trials) %>% 
  mutate(trial = factor(trial,
                        levels = trials,
                        labels = str_c("trial ", labels)))
  

func_load_image = function(situation){
  readPNG(str_c("../../figures/stimuli/trial_", situation, ".png"))
}

# linking images and clips
df.trials = df.plot %>% 
  distinct(trial) %>% 
  arrange(trial) %>% 
  mutate(number = trials,
         grob = map(.x = number, .f = ~ func_load_image(situation = .x)),
         label = str_c("trial ", number))

df.text = df.plot %>%
  distinct(trial) %>% 
  arrange(trial) %>% 
  mutate(x = -Inf,
         y = 1.2,
         label = 1:n())

# plotting
ggplot(data = df.plot,
       mapping = aes(x = item,
                     y = value)) +
  geom_col(data = df.plot %>% 
             filter(index == "people"),
           mapping = aes(fill = question_type,
                         linetype = question_type),
           alpha = 0.5,
           color = "black") +
  geom_linerange(mapping = aes(ymin = low,
                               ymax = high),
                 size = 1.5) + 
  geom_point(data = df.plot %>% 
               filter(index %in% c("rationality", "pivotality", 
                                   "mixture")) %>% 
               mutate(index = factor(index, levels = c("rationality", "pivotality", 
                                   "mixture"))),
             mapping = aes(shape = index,
                           fill = index),
             position = position_dodge(width = 0.5),
             color = "black",
             stroke = 1,
             size = 4) +
  geom_custom(data = df.trials,
              mapping = aes(data = grob,
                            x = -Inf,
                            y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.05,
                                                hjust = 0)) + 
  geom_text(data = df.text,
            mapping = aes(x = x,
                          y = y,
                          label = label),
            hjust = -0.5,
            size = 20,
            color = "white") +
  facet_wrap(~trial,
             labeller = labeller(label),
             scales = "free_x",
             nrow = 2) +
  labs(linetype = "question",
       shape = "model") + 
  coord_cartesian(clip = "off",
                  ylim = c(0, 1)) + 
  scale_y_continuous(breaks = seq(0, 1, 0.25),
                     labels = seq(0, 100, 25),
                     expand = expansion(add = c(0, 0))) +
  scale_shape_manual(values = c(21, 23, 24)) + 
  scale_linetype_manual(values = rep(1, 3)) + 
  scale_fill_manual(breaks = c("choice", "strength", "trees",
                               "rationality", "pivotality", "mixture"),
                    values = c(choice = "gray80", 
                               strength = "gray40", 
                               trees = "gray0", 
                               rationality = "blue", 
                               pivotality = "red", 
                               mixture = "purple")) + 
  theme(legend.position = "bottom",
        legend.title = element_text(size=30),
        legend.text = element_text(size=25),
        panel.grid.major.y = element_line(),
        axis.text.y = element_text(size = 25),
        axis.text.x = element_text(size = 25),
        axis.title = element_blank(),
        strip.background = element_blank(),
        strip.text = element_blank(),
        panel.background = element_rect(fill = NA, color = "black"),
        panel.spacing.x = unit(0.5, "cm"),
        panel.spacing.y = unit(7, "cm"),
        plot.margin = margin(t = 7, l = 0.2, r = 0.2, b = 0.1, unit = "cm")) +
  guides(fill = F,
         linetype = guide_legend(override.aes = list(fill = c("gray80",
                                                              "gray40",
                                                              "gray0"))),
         shape = guide_legend(override.aes = list(fill = c("blue", "red", "purple"), 
                                                  size = 8))) 
  
# ggsave("../../figures/plots/selection.pdf",
#        width = 20,
#        height = 13)

5 Session info

R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] forcats_0.5.1     stringr_1.4.0     dplyr_1.0.4       purrr_0.3.4      
 [5] readr_1.4.0       tidyr_1.1.2       tibble_3.0.6      tidyverse_1.3.0  
 [9] rsample_0.0.9     modelr_0.1.8      lubridate_1.7.9.2 kableExtra_1.3.1 
[13] ggtext_0.1.1      patchwork_1.1.1   egg_0.4.5         ggplot2_3.3.3    
[17] gridExtra_2.3     png_0.1-7         xtable_1.8-4      DT_0.17          
[21] corrr_0.4.3       tidyjson_0.3.1    RSQLite_2.2.3     janitor_2.1.0    
[25] knitr_1.31       

loaded via a namespace (and not attached):
 [1] nlme_3.1-151      fs_1.5.0          bit64_4.0.5       webshot_0.5.2    
 [5] httr_1.4.2        tools_4.0.3       backports_1.2.1   R6_2.5.0         
 [9] DBI_1.1.1         mgcv_1.8-33       colorspace_2.0-0  withr_2.4.1      
[13] tidyselect_1.1.0  bit_4.0.4         compiler_4.0.3    cli_2.3.0        
[17] rvest_0.3.6       xml2_1.3.2        bookdown_0.21     scales_1.1.1     
[21] digest_0.6.27     rmarkdown_2.6     pkgconfig_2.0.3   htmltools_0.5.1.1
[25] parallelly_1.23.0 dbplyr_2.0.0      fastmap_1.1.0     highr_0.8        
[29] htmlwidgets_1.5.3 rlang_0.4.10      readxl_1.3.1      rstudioapi_0.13  
[33] farver_2.1.0      generics_0.1.0    jsonlite_1.7.2    crosstalk_1.1.1  
[37] magrittr_2.0.1    Matrix_1.3-2      Rcpp_1.0.6        munsell_0.5.0    
[41] lifecycle_1.0.0   furrr_0.2.2       stringi_1.5.3     yaml_2.2.1       
[45] snakecase_0.11.0  blob_1.2.1        parallel_4.0.3    listenv_0.8.0    
[49] crayon_1.4.1      lattice_0.20-41   haven_2.3.1       splines_4.0.3    
[53] gridtext_0.1.4    hms_1.0.0         pillar_1.4.7      markdown_1.1     
[57] codetools_0.2-18  reprex_1.0.0      glue_1.4.2        evaluate_0.14    
[61] vctrs_0.3.6       cellranger_1.1.0  gtable_0.3.0      future_1.21.0    
[65] assertthat_0.2.1  cachem_1.0.1      xfun_0.21         broom_0.7.3      
[69] viridisLite_0.3.0 memoise_2.0.0     globals_0.14.0    ellipsis_0.3.1